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We study the phase diagram of mass- and spin-imbalanced unitary Fermi gases, in search for 
the emergence of spatially inhomogeneous phases. To account for fluctuation effects beyond the 
mean-field approximation, we employ renormalization group techniques. We thus obtain estimates 
for critical values of the temperature, mass and spin imbalance, above which the system is in the 
normal phase. In the unpolarized, equal-mass limit, our result for the critical temperature is in 
accordance with state-of-the-art Monte Carlo calculations. In addition, we estimate the location 
of regions in the phase diagram where inhomogeneous phases are likely to exist. We show that 
an intriguing relation exists between the general structure of the many-body phase diagram and 
the binding energies of the underlying two-body bound-state problem, which further supports our 
findings. Our results suggest that inhomogeneous condensates form for mass imbalances m 4 ./m^ > 3. 

The extent of the inhomogeneous phase in parameter space increases with increasing mass imbalance. 


I. INTRODUCTION 

The past 15 years have witnessed tremendous advances 
in the experimental control and exploration of ultracold 
atomic Fermi gases. Since the first realization of a Bose- 
Einstein-condensate of paired fermions [T], experimental 
techniques have been further developed and now allow 
for detailed studies of many-body phenomena, controlled 
variations in temperature and polarization [5], studies of 
Bose-Fermi mixtures [3], optical lattices [1], as well as 
precise determinations of the equation of state [5], see 
Refs. [6H8] for reviews. It has thus become possible to 
test theoretical descriptions of long-known effects such as 
Bardeen-Cooper-Schrieffer (BCS) superfluidity or Bose- 
Einstein-condensation (BEC) with high precision. More¬ 
over, many new experimental studies of many-body phe¬ 
nomena with mixtures of a variety of different fermion 
species (such as ®Li, ^®^Dy, ^®^Dy, and ^®^Er) are 
within reach in the near future (see e.g. Ref. [QlUljl. 
giving us an unprecedented opportunity to better our 
understanding of mass imbalances in strongly coupled 
Fermi gases and push our understanding of more ex¬ 
otic phenomena such as the emergence of inhomogeneous 
phases [HUS] to a whole new level. 

In experiments, the particle density n and s-wave scat¬ 
tering length tts are control parameters. In a sufficiently 
dilute gas, they represent the only scales of the systems, 
since the effective range of the interaction can safely 
be neglected. Experimentally, Qs can be tuned at will by 
means of so-called magnetic Feshbach resonances. This 
opens up the possibility to explore many-body phenom¬ 
ena over a wide range of interaction strengths. Of par¬ 
ticular interest is the strongly coupled “unitary” limit, 
where Og —>■ oo. In this regime, a non-perturbative treat¬ 
ment is inevitable due to the absence of a small expansion 
parameter |14j . 

Systems with equal populations and particle masses 
for the different fermion species, i.e. spin- and mass- 
halanced unitary Fermi gases, are well under control by 
now from both an experimental and theoretical point 


of view. For example, lattice Monte Carlo studies of 
the equation of state and the critical temperature have 
reached a quantitative level m and show good agree¬ 
ment with experimental data HU. For spin- and mass- 
imbalanced Fermi gases, which are at the heart of this 
work, less is known beyond the mean-field approxima¬ 
tion, although great efforts have been made in recent 
years to study mass-imbalanced (see, e.g., Refs. [IM]) 
as well as spin-imbalanced (see, e.g. Refs. [33H3T] 1 uni¬ 
tary Fermi gases (see, e.g.. Refs. [32l [33] for reviews). 
The difficulties encountered in studies of imbalanced sys¬ 
tems beyond the mean-field limit are many. For example, 
ab initio Monte Carlo studies are severely hampered by a 
so-called sign problem if spin and/or mass imbalances are 
introduced [53J|35]. Techniques to surmount this problem 
have been developed [35], but have so far focused on the 
zero-temperature equation of state of mass-imbalanced 
unitary Fermi gases m, and their use is very recent. 

The reasons for the interest in imbalanced systems are 
manifold as well. For example, the Fermi surfaces of the 
different species are mismatched in this case, possibly 
giving rise to exotic phenomena such as inhomogeneous 
phases of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) 
type [T2|[T3] or Sarma [38| phases. For one-dimensional 
spin-imbalanced Fermi gases, where inhomogeneous pair¬ 
ing is expected to exist for a wide range in parameter 
space [33], observation of an FFLO phase has indeed 
been claimed [40] . In three dimensions, however, the 
inhomogeneous phase is expected to occupy only a thin 
layer of parameter space between the homogeneous su¬ 
perfluid and the normal fluid in parameter space, if at 
all [4lti43] . This renders the experimental detection of 
such phases quite challenging. The utilization of mass- 
imbalanced mixtures is expected to alleviate this situ¬ 
ation somewhat due to the larger parameter space for 
inhomogeneous pairing in this case |44j . 

Most of the studies, especially of FFLO phases, have 
so far relied on the mean-field approximation. However, 
even in three dimensions, these studies yield at best qual¬ 
itative insights into the phase structure of the system. In 
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fact, even in the balanced case it is known that the crit¬ 
ical temperature (measured in units of the chemical po¬ 
tential) is overestimated by a factor ^ 1.6 in mean-field 
studies. This can be traced back to the omission of order- 
parameter fluctuations. However, the fate of FFLO-type 
phases upon inclusion of such fluctuation effects remains 
largely an open question. 

In this work, we analyze the phase structure of a uni¬ 
tary Fermi gas with spin and mass imbalance at finite 
temperature. For this purpose, we employ functional 
renormalization group (fRG) techniques, which allow us 
to include order-parameter fluctuations. Results from 
previous studies for the superfluid critical temperature 
of the mass- and spin-imbalanced case are in good agree¬ 
ment with experimental data and ab initio MC stud¬ 
ies [3T]. Here, we extend the theoretical framework de¬ 
veloped in earlier RG works [STJ I4M48] in order to inves¬ 
tigate mass- and spin-imbalanced systems. Although our 
setup does not yet allow us to explicitly resolve inhomoge¬ 
neous phases, strong hints of their existence beyond the 
mean-field approximation can already be detected and 
we discuss them here. 

The rest of the paper is organized as follows. In 
Sec. we introduce the microscopic model and define 
the scales and dimensionless parameters. Since the iden¬ 
tification of inhomogeneous phases is challenging, we be¬ 
gin our discussion of the phase structure of mass- and 
spin-imbalanced Fermi gases by revisiting the two-body 
bound-state problem in Sec. |HI| In a study of one¬ 
dimensional gases [33], it was found that the two-body 
problem in the presence of Fermi spheres provides useful 
information about the many-body problem. In fact, it 
was then found that the phase structure is in approxi¬ 
mate quantitative agreement with the associated many- 
body study in the mean-field approximation. In Sec. |IV| 
we present the mean-field phase diagram. Gorrections 
beyond the mean-field approximation are discussed in 
Sect. jVj and the phase diagram resulting from those cor¬ 
rections is shown. 


s-wave scattering length Og: 
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where A is the ultraviolet cutoff and Creg is a constant 
that depends on the regularization scheme. 

It is convenient to trade in the two fermion masses 
for an imbalance parameter m using the following defini¬ 
tions: 


m+ = -—^ , m_ 
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In this work we set to+ = 1 which corresponds to the 
choice = 1/2 in the mass-balanced case. Note that 
the mass imbalance parameter to is thus normalized such 
that 0 < TO < 1. The chemical potentials of the fermion 
species can be expressed in terms of an average chemical 
potential ^ and the (normalized) spin imbalance param¬ 
eter or so-called Zeeman field h via 

( 4 ) 

It is also convenient to define a dimensionless tempera¬ 
ture parameter 
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Finally, we choose units such that /i = fee = 1- 

The action S'p defined above features a global U{1) 
symmetry associated with particle number conserva¬ 
tion. This symmetry is spontaneously broken by a non¬ 
vanishing field expectation value (V’tf/’;) ^ 0 in the su¬ 
perfluid phase. The main goal of the present work is 
to identify the region of parameter space where 
as a function of {h, rh, T), becomes nonzero and possibly 
position-dependent, indicating the breakdown of transla¬ 
tional invariance. 


II. MICROSCOPIC MODEL 

Microscopically, the two-component spin- and mass- 
imbalanced Fermi gas in the vicinity of a broad Feshbach 
resonance is described by the following action: 

.t?:; V 2 to. ; 

Here, g = f dr f (Px. The two fermion species are rep¬ 
resented by the Grassmann-valued fields 'ipa- that depend 
on spatial coordinates x and compact Euclidean time r. 
A contact interaction of strength g couples both species. 
Its strength can be tuned through its dependence on the 


III. TWO-BODY BOUND STATES 

It is particularly challenging to identify those regimes 
in the phase diagram where the C/(l) symmetry and 
translational invariance are broken simultaneously. For 
example, it may be reasonable to assume that the order 
parameter {%p^ip^){x) is a periodic function. However, 
neither its precise functional form nor the length scale 
associated with its period (i.e. the characteristic mo¬ 
mentum Q associated with the inhomogeneity in Fourier 
space) are known a priori. Therefore, it is worthwhile to 
perform preparatory analyses to help identify domains in 
parameter space where U{1) symmetry breaking may ap¬ 
pear in the full many-body problem and, in particular, 
where such breaking is accompanied by a spontaneous 
breakdown of translation invariance. 

The physical interpretation of a finite order parameter 
7 ^ 0 is a condensate of paired fermions. Naturally, 
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those (bosonic) pairs have to be formed in the first place 
in order to condense. It thus seems reasonable to assume 
that the characteristics of the energetically most favor¬ 
able pairing state (for a given parameter set (h,fh,T)) 
will strongly influence the properties of a (potential) con¬ 
densate in the many-body problem. A quantum mechani¬ 
cal bound-state calculation in Ref. [51] has indeed shown 
good qualitative and even semi-quantitative agreement 
with the many-body mean-field analysis in one dimen¬ 
sion. Here, we perform a similar calculation for the three- 
dimensional case. 

Following Ref. |39|, the Schrodinger equation for the 
wave function 'I' of two distinct fermions in the presence 
of their respective Fermi surfaces is given by 



X! - xi) -f Eb 




cr=t,t 

( 6 ) 

where Eb = F' and the delta-shaped potential 

is the two-body equivalent of the contact interaction in 
Eq.0. Thus, the relation between the coupling strength 
g and the s-wave scattering length is given by Eq. ([^. 
The kinetic energy is measured with respect to the Fermi 
surfaces: 6^(9^,^) = | - - ep.trl, where eF,,^ 

corresponds to the Fermi energies of the two species, re¬ 
spectively. Note that the general setup is well known and 
has been previously used to determine the properties of a 
single Cooper pair in the context of balanced BCS theory, 
see e.g. Refs. [ 511155 ] . 

In momentum space, Eq. <§ can be recast into a 
(renormalized) integral equation for the binding en¬ 
ergy Eb- 
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(l+pj 

+ 

(l-pj 

+ Eb{\P\) 



where f^=f d^p/(27r)^. The vectors p and P denote 
the relative and center-of-mass momenta, respectively. 
Since the energy of the two-particle state 'k depends on 
the magnitude of the center-of-mass momentum P, the 
ground-state solution of Eq. 0 determines whether a 
bound state exists for a specific set of parameters and, 
crucially for this work, provides information about the 
energetically favored center-of-mass momentum. 

In Fig.[^we show the results for the ground-state bind¬ 
ing energies. We find that for small m, and around the 
line of equal Fermi momenta fcp.t = ^F.i at h = ?n, 
the bound-state formation with zero center-of-mass mo¬ 
mentum is favored (light-gray shaded area). Due to 
the increasing mismatch of Fermi surfaces away from 
the h = m line, the (dimensionless) binding energy 
Eb = Eb/p decreases monotonically until it reaches 
zero for high mass imbalances and negative spin imbal¬ 
ances. If, on the other hand, fh > 0.46, pairing with 
finite P = P /is favored for sufficiently small spin im¬ 
balances (dark/red-shaded area). Note that this seems 


Figure 1. Dimensionless two-particle binding energies Eb = 
Eb/ p as obtained from Eq. 0 - Domains in parameter space 
where bound states are found {Eb > 0) are (gray) shaded. 
For red/dark shading, bound-state formation with a finite 
center-of-mass momentum P = |P|/ydI is favored. The 
line h = fh corresponds to the case of equal Fermi momenta 
for the two species. The two-body states found along this line 
are found to be the overall most deeply bound states. 


to stabilize the binding energy away from equal Fermi 
momenta for very large mass imbalances, manifested by 
a slight back-bending of the Eb isolines. Similar behav¬ 
ior was observed in the one-dimensional case [3Hj . where 
it was found to have strong influence on the structure of 
the many-body phase diagram. 

Loosely speaking, a condensate of pairs with finite 
center-of-mass momentum would break translational in¬ 
variance. Our observation of the existence of a region of 
parameter space associated with two-body bound states 
with a finite center-of-mass momentum suggest that a 
region governed by an inhomogeneous ground state may 
exist in the many-body phase diagram. 

We close this section with a word of caution as to 
the relevance of our two-body analysis for the actual 
many-body problem: The existence of bound states in 
the two-body problem in the presence of Fermi surfaces 
does not necessarily entail spontaneous symmetry break¬ 
ing in the associated many-body problem. The latter re¬ 
quires, additionally, Bose-Einstein condensation of said 
bound states. Furthermore, the consideration of inert 
Fermi surfaces in our two-body study is questionable in 
the strongly coupled unitary regimeFlFor example, the 
so-called Fermi-polaron problem m [24l - [30] . which con¬ 
stitutes the limiting case of a single spin-up impurity in 
a sea of indistinguishable spin-down fermions, is known 


^ Strictly speaking, the assumption of inert Fermi surfaces is only 
justified if the Fermi momenta/energies of the non-interacting 
system entering our two-body study are of the order of the Fermi 
momenta/energies attributed to the fully interacting many-body 
problem. 
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to have a negative chemical potential. This implies that 
energy is gained when a spin-up impurity is added to a 
system of non-interacting spin-down fermions. As dis¬ 
cussed in Ref. [53] , this energy gain determines the value 
of the Zeeman field h above which a mixed-phase may 
emerge in experimental studies. Moreover, this value 
of h can be viewed as a strict lower bound for the critical 
value of h above which the ground state of the system 
is superfluid. Note that the interaction of the impurity 
with all the spin-down fermions is taken into account in 
these studies. Our study of the two-body problem in the 
presence of Fermi spheres cannot reproduce the results 
of the above-mentioned Fermi-polaron problem, as we 
only allow for an interaction between one spin-up and one 
spin-down fermionj^ Our approach should therefore be 
viewed as complementary to the Fermi-polaron problem. 
In fact, our analysis gives direct access to the center-of- 
mass momentum of the bound state and therefore enables 
us to estimate the regions in the many-body phase dia¬ 
gram where inhomogeneous phases are likely to exist. In 
any case, all results so far are by construction valid only 
at strictly zero temperature. Therefore, a proper many- 
body treatment is mandatory in order to obtain the ac¬ 
tual phase diagram. However, as we will see below, our 
predictions resulting from Eq. Q for the position of in¬ 
homogeneous phases turn out to be astonishingly good, 
which emphasizes the importance of few-body physics for 
our understanding of complex many-body phenomena. 


IV. MEAN-FIELD ANALYSIS 


A. Formalism and Fulde-Ferrell Ansatz 


For this work, the properties of the order parame¬ 
ter are essential. We therefore formulate the 

problem of spontaneous symmetry breaking in terms 
of the associated order-parameter field. This can be 
achieved by means of a judiciously chosen Hubbard- 
Stratonovich transformation EH, upon which we obtain 
the following microscopic action 


>S'b[{V'<t},V5] 


a —, 4 " 






■ ( 8 ) 


This action is equivalent to the purely fermionic ac¬ 
tion Sp in Eq. Q and is the basis for all the calcu¬ 
lations in this work. The boson field (p ^ 
be viewed as mediating the fermionic self-interaction 
~ Note that no kinetic term 


^ Note also that the non-interacting Fermi spheres enter our com¬ 
putation and that the associated Fermi momentum kp -f- (fcp 
becomes complex for h < —1 {h > 1). 


for the order-parameter field appears in the classical ac¬ 
tion ‘5 'b- Such a term is generated dynamically when 
fermion fluctuations are integrated out, due to the pres¬ 
ence of the Yukawa-type interaction, as we will see in 
more detail below. 

From a field theoretical point of view, the field (p 
is nothing but an auxiliary field, introduced by the 
Hubbard-Stratonovich transformation to facilitate the 
computations by removing the quartic fermion term in 
favor of a Yukawa-type interaction. From a phenomeno¬ 
logical point of view, however, (p may be interpreted 
as a collective state of two fermions, corresponding to 
the closed channel of the Feshbach resonance, see, e.g.. 
Refs. ggEg. 

In our partially bosonized formulation, the order pa¬ 
rameter for U{1) symmetry breaking can now be iden¬ 
tified with ipo ~ In the ground state, these 

field expectation values can be obtained by minimizing 
the quantum effective action F ^ — In 2^ with respect to 
ip ^ where 

Z = J = J -Dpdet^lip] (9) 

is the partition function in the path-integral representa¬ 
tion. Note that the determinant on the right-hand side 
can in principle be used to define a purely bosonic ac¬ 
tion S'pB [<^] = — In det^ [ip ], which is the starting point 
for lattice MC calculations [TSl IMl 155] (see Ref. [33] for 
a review). While the effective action F shares the s;^- 
metries of the microscopic action S-q by construction!^ it 
includes the effects of all thermal and quantum fluctua¬ 
tions. Thus, it is composed of renormalized fields and 
couplings which are in general momentum dependent. 
Since the fermion determinant det^ involves a generally 
complicated dependence on p, an exact computation of 
F is highly non-trivial. Therefore, systematic approxima¬ 
tion schemes are needed to gain insight into the physical 
content of the theory. 

In this sense, the widely used mean-field approxima¬ 
tion can be considered as a lowest-order approximation 
to the effective action. It is obtained by shifting the 
field Cp ^ (p -\- 6ip, where 5(p now represents the fluctua¬ 
tion field, and performing a saddle-point approximation 
of the path integral about (p. This renders the analytic 
computation of the fermion determinant more feasible. 
However, it is still non-trivial to carry out if one allows 
for a general space-dependent “background” field (p{x), 
which is needed to enable the detection of inhomogeneous 
phases. In the following, we employ an analytically fea¬ 
sible Fulde-Ferrell ansatz (FF) [12] . 

ip{x) = . (10) 


Since ip is chosen to be a constant amplitude, we are left 
with the standard mean-field ansatz in the limit of van¬ 
ishing momentum Q. Note that the ansatz (10) may be 


® The path integral measure respects the symmetries of the theory. 
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regarded as the first term of a Fourier expansion of a 
more general (p{x). It is hence considered to be a par¬ 
ticularly good approximation to the full solution in the 
vicinity of a second order phase transition to the normal 
phase, where higher order contributions are expected to 
be small (see e.g. Refs. [Ml [54U56] i. 

Using the ansatz (10), the order-parameter potential U 
becomes 

U{^tQ) = j — E^ ' 


2g2 


a=±l 


Tlog(l 




where 


= V ( 9 ^ +Q^ - 2mq- Q - /r)^ -|- , 


( 11 ) 


( 12 ) 


and Q = \Q\ as well as q= The quantity A = 
evaluated at the (global) minimum (pQ of U is nothing 
but the fermion gap. Note that the minimum of U in 
Q-direction is not necessarily an extremum. A global 
(numerical) minimization of the potential is therefore re¬ 
quired to find the ground state of the theory. 


B. Results 


In Fig. we show the phase diagram as obtained from 
a direct minimization of U{A,Q) in Eq. (111. Here, we 
mainly discuss aspects related to the emergence of inho¬ 
mogeneous phases. For detailed discussions on the extent 
and the properties of the homogeneous phases at zero and 
finite temperature, we refer the reader to earlier work, 
e.g. Refs. 

We begin by briefly discussing some characteristic fea¬ 
tures of the phase diagram. The region of homogeneous 
symmetry breaking is roughly centered around the line 
of equal Fermi momenta, h = rh, as su gges ted by the 
analysis of the two-body problem of Sec. O At T = 0, 
the transition from the homogeneous superfluid to the 
normal phase is always of first order for the parameter 
space considered in the present work. At sufficiently high 
temperatures, on the other hand, the transition from the 
superfluid to the normal phase changes into a second- 
order transition. The surfaces in parameter space associ¬ 
ated with these two types of transitions meet at a critical 
line denoted by T),p. Note that we limit ourselves to spin 


imbalances h G [—1,1]. Contrary to the case of the two- 
body analysis in Sec. |III[ this constraint is not imposed 
by the method itself. In fact, the domain of high mass 
imbalance, m > 0.8 for /i > 1 is of interest for an investi¬ 
gation of the physics of Sarma phases |35]. In particular, 
it was found in mean-field calculations that a second or¬ 
der transition occurs in this regime even at T = 0 |18j[^ 


By employing the functional RG scheme discussed in Sec. ly^ 



Figure 2. Mean-field phase diagram in the plane spanned 
by h and m with T-isolines. The (light) gray-shaded region 
corresponds to a homogeneous condensate, i.e. Aq > 0 and 
Qo = 0 (SF); the dark/red-shaded area depicts a Fulde-Ferrell 
regime, i.e. Aq > 0 and Qo > 0 (FF). The (blue) dotted lines 
correspond to lines of multicritical points TJ-p, where either 
a surface associated with first-order transitions from the ho¬ 
mogeneous to the inhomogeneous phase meets a surface as¬ 
sociated with second-order transitions from the homogeneous 
to the normal fluid (NF) phase as well as a surface from the 
inhomogeneous to the normal fluid phase (NF), or where a 
surface associated with second-order transitions from the ho¬ 
mogeneous to the NF phase meets a surface associated with 
first-order transitions from the homogeneous to the NF phase. 


We leave a detailed discussion of Sarma phases for future 
work. 

We now turn to a discussion of the inhomogeneous 
phase depicted by the dark/red-shaded area in Fig. 
Studies similar to ours have been performed for specific 
values for m (e.g., Li-K mixture |23j ) as well as arbitrary 
values of rh employing a T-matrix approach [44j . While 
our mean-field results agree well with those of Ref. [53] , 
we do not find an inhomogeneous phase for mass im¬ 
balances down to m = 0 as in Ref. [44]. There is, in 
fact, disagreement in the literature on this point (see e.g.. 
Refs. HmH]), which suggests that the appearance of an 
inhomogeneous phase is very sensitive to the details of 
the approximation scheme as well as the specific ansatz 
for the inhomogeneity. This strengthens our motivation 
to consider alternative sources of information on pairing 
behavior such as our few-body analysis of Sec .|III| and RG 
methods to account fluctuation effects (Sec. |V| below). 

Fig.|3| displays an enlarged view of the region with 
inhomogeneous superfluidity, focusing on the domain of 
parameter space governed by inhomogeneous pairing at 


we do indeed find strong hints that a zero-temperature Sarma 
phase in the regime of large m and h > 1 exists even beyond the 
mean-field approximation. For a discussion of the existence of 
Sarma phases in the phase diagram of spin-imbalanced but mass- 
balanced Fermi gases beyond the mean-field limit see Ref. ESI- 
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Figure 3. Zoom of Fig. showing the regime in parame¬ 
ter space where inhomogeneous phases (FF) are found in our 
mean-held study. The multicritical (dotted) and the second- 
order (dashed and dot-dashed) lines are dehned as in Fig. 
and included here to guide the eye. The black/solid line rep¬ 
resents the transition line between the homogeneous (SF) and 
inhomogeneous or normal phases (NF) at T = 0, respectively. 
Finally, the (red) dot-dot-dashed line is the lower bound for 
pairing with P > 0 as obtained from our two-body analysis, 
see Fig. 


f = 0. As evident from the figure, the inhomogeneous 
phase occupies only a thin shell in parameter space close 
the phase with a homogeneous condensate. It is sepa¬ 
rated from the normal phase by a second-order transition 
and from the homogeneous phase by a first-order tran¬ 
sition]^ Note that, for finite negative scattering lengths 
Og, a similar behavior has been found previously also for 
the mass-balanced case [32] • 

Within numerical errors, the position of the first-order 
transition line from the homogeneous to the inhomoge¬ 
neous phase coincides almost perfectly with the transi¬ 
tion line from the homogeneous superfluid to the nor¬ 
mal phase, which is found if one does not allow for a 
space-dependent order parameter]^ The observed first- 
order nature of the homogeneous-inhomogeneous transi¬ 
tion may be an artifact of the FF-ansatz in our study. 


® Both the phase with a homogeneous condensate and the one 
with an inhomogeneous condensate are associated with sponta¬ 
neous C/(l) symmetry breaking. The boundary between these 
two phases is solely associated with translation symmetry break¬ 
ing where the value of Q = |Q| of the condensate acts as an 
order parameter. On the other hand, C/(l) symmetry is restored 
in the normal phase. The boundary between the normal phase 
and the phase with a (in)homogeneous condensate is therefore 
associated with 17(1) symmetry breaking where Cpo plays the role 
of the order parameter. 

® For fixed fh, we find that the extent of the homogeneous phase 
which has been superseded by an inhomogeneous phase is as 
narrow as 5h ^ 0.006 for the largest considered m = 0.95 and 
decreases rapidly towards the zero-temperature endpoint of the 
critical line at CP® ~ (0.09,0.65). 


In fact, studies of analytically solvable relativistic mod¬ 
els in ID found that the corresponding transition line is 
of second order rather than first order [SU |SS] . The po¬ 
sition of the associated transition line was found to be 
reasonably well described when the space dependence of 
the condensate is approximated by the first term of the 
Fourier expansion [56] . as done by the FF ansatz. How¬ 
ever, the order of the transition is spuriously predicted 
to be of first oder0 

The very small extent of the inhomogeneous phase in 
h-direction for fixed fh close to the zero-temperature crit¬ 
ical endpoint CP° renders the precise determination of its 
coordinates very difficult from a numerical point of view. 
Following our line of arguments from Sec. |III[ the quick 
disappearance of the inhomogeneous phase close to this 
point is not unexpected: The two-body binding energy 
away from its maximum at h = m is “stabilized” for finite 
pair-momentum only at very large mass imbalances, as 
indicated by the characteristic back-bending of the lines 
of constant binding energy in Fig. Since deeply bound 
pairs should be less sensitive to (quantum) fluctuations 
that tend to destroy pairing correlations, the formation of 
a condensate in the many-body problem may be expected 
to be favored in this regime as well. Thus, the inhomo¬ 
geneous phase is expected to widen towards smaller h 
only at large fh, whereas it is expected to narrow down 
when fh is decreased. This is indeed what we find in 
Fig.jl 

In the same figure, we also show the lower bound (dot- 
dot-dashed line) for pair formation with a finite center- 
of-mass momentum P > 0, see also Fig. We find that 
there is an obvious discrepancy between this line (which 
we argued above is important for the occurrence of in¬ 
homogeneous condensates) and the actual many-body 
phase boundary associated with a transition from a ho¬ 
mogeneous to an inhomogeneous phase. 

The critical endpoint CP° of the inhomogeneous phase 
appears to coincide quite well with the intersection point 
of the lower bound for finite-momentum pair formation 
(dot-dot-dashed line) and the phase boundary of the 
phase with a homogeneous condensate (black solid line). 
A high precision determination of the position of CP° 
is beyond our reach at the moment due to the numeri¬ 
cal complications mentioned above; we therefore refrain 
from drawing rigorous conclusions about this stunning 
coincidence. For larger values of fh and h in the domain 
between the (black) solid and (red) dot-dot-dashed line, 
we find that a homogeneous condensate is energetically 
favored over an inhomogeneous condensate (see Fig. |^, 
even though bound-state formation with finite center-of- 
mass momentum should still be preferred. However, as 
discussed above, the simple FF ansatz (10) for (f{x) is not 
expected to yield precise results away from the transition 


^ Note that, in addition to a first-order transition measured by 
the parameter Q, we observe a discontinuity in the (7(1) order- 
parameter ifio as well. 
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to the normal phase. While our two-body calculation is 
not hampered by assumptions of this kind, it is limited 
by the assumption of inert Fermi surfaces. In contrast, 
our mean-field study of the full many-body problem takes 
into account interactions between all spin-up and spin- 
down fermions. We shall come back to this issue below 
when we discuss the results of our RG study. 

In summary, the need to investigate the dynamics 
of the many-body problem beyond the mean-field limit 
seems inevitable. Furthermore, the above analysis indi¬ 
cates that one must include more sophisticated ansatze 
for the condensate (p{x) to gain reliable insights into the 
problem of inhomogeneous superfluidity. Our fRG anal¬ 
ysis, as detailed next, aims to provide a systematic way 
to approach this goal. 


V. BEYOND THE MEAN-FIELD 
APPROXIMATION 

A. Functional Renormalization Group: Formalism 


with renormalized quantities 




_ 1 _ 
u — A 2 u 
fL^ — il^ip • 


(16) 


This ansatz is an extension of a class of ansatze which 
has been successfully applied to spin-balanced (see, e.g.. 
Refs. 122 122 ]) spin-imbalanced (see, e.g.. 

Refs. |30ll3ll[59l[22]) systems across the whole BEG-BCS 
crossover. 

In addition to an ansatz for F^, the solution of the RG 


equation (13) requires an initial condition F^^a- In our 


case this implies choosing initial conditions for the so- 
called wavefunction renormalizations and the 

Yukawa coupling h^, and the effective order-parameter 
potential Uk- The latter depends only on the C/(l) in¬ 
variant p = If*If. At the UV scale A, it is fixed by the 
two-body problem in the vacuum and is therefore simply 
given by [45] 


Ua{p) = 


(17) 


where 


To understand spontaneous symmetry breaking in 
quantum systems, it is essential to account for quantum 
fluctuations of the order-parameter. In order to accom¬ 
plish this in a systematic way, we employ a fRG approach. 
More specifically, the Wetterich equation [20] 


SfeFfe = ^STr 


dkRk 


(13) 


will be central to our analysis. Here, k is the RG scale in¬ 
troduced by the regulator function Rk- The latter spec¬ 
ifies the details of the Wilsonian momentum-shell inte¬ 
gration. In particular, it provides for a suitable infrared 
(IR) and ultraviolet (UV) regularization, see App. [^for 
details. The Wetterich equation determines the change 
of the effective average action F^ under a change of the 
scale k and therefore allows to interpolate between the 
microscopic action S' at a given UV cutoff scale A and the 
full quantum effective action F = F k^o in the long-range 
(i.e. IR) limit: 


S = lim Ffe, F = lim Ffc. (14) 

fc—>-A—>-oo fc—>-0 


For reviews and introductions to this approach with re¬ 
spect to an application to ultracold Fermi gases, see 
Refs. 1^. 

While Eq. (13) is exact, in general one must consider 
an ansatz for F^ (i.e. a truncation of the full F^) in order 
to arrive at a solution. Here, we choose 


^k[{M,f] = [ \ 

j T.X ^ _ .A. I ^ 


dr- 


f 


( 


Zp^k 


A. 


(p,k 


dr- 


2m„ 


l-fh^ 


V'o 


0 


‘P 


(15) 


Uk{p) - hp[f*'tp^'tl)^ - f^ip*^ 


_ I A 


(18) 


and 


hp^A = VGtt^A . 


(19) 


Here, the initial value of the Yukawa coupling has been 
chosen such that the system is set to be close to a broad 
Feshbach resonance [ 22 ]. The parameter ~ {B — Bq) 
measures the detuning of the system with respect to the 
resonancej^ Since we are only interested in the unitary 
limit (i.e. the system at resonance), we set i^a = 0 for 
the remainder of this work. Higher order terms ^ p" in 
the potential Uk are generated during the RG flow due to 
quantum fluctuations and are taken into account in our 
analysis. 

At this point we are left with the determination of 
the initial conditions for the wavefunction renormaliza¬ 
tions. The momentum and frequency dependence of 
the boson-field propagator effectively allows to resolve at 
least part of the momentum and frequency dependence 
of the fermionic self-interactions (see, e.g.. Ref. 161) for 
a more general introduction). Per our ansatz (15), the 
inverse boson propagator for the bare/unrenormalized 
field f in momentum space i^ 




(^ZtpkQO 


l — m?^ ^2 , 

2 dfdf* 


), ( 20 ) 

I' 


® Here, B denotes the external magnetic field. 

^ For convenience, we do not show contributions to the propagator 
stemming from the insertion of the regulator function into the 
path integral. 



















where ipo denotes the /c-dependent value of the position 
of the ground-state of the potential Uk- The wavefunc- 
tion renormalizations and are assumed to be 
independent of and q. Here, we only take into ac¬ 
count the scale dependence of the wavefunction renor¬ 
malizations, which is minimally required to detect the 
emergence of an inhomogeneous ground state within our 
present setup (see our discussion below). The structure 
of this propagator can be understood as arising from a 
derivative expansion of the effective action that has been 
truncated at the lowest nontrivial order. 

The RG flow of the wavefunction renormalizations Z^p^k 
and is conveniently parameterized with the aid of 
the associated anomalous dimensions rjz^k and rjA,k, re¬ 
spectively. For example. 


'nA,k = -kdk\^Ap,^k- ( 21 ) 


In the present work, we keep the ratio Z^^k/A^^k = 1 
fixed in the RG flow and compute only the flow of Ap,^k, 
see also Ref. [3T] . The initial condition for A^ k is given 
by 

lim Ap^k = 0 . (22) 

k—^A—^oc 


From a phenomenological point of view, this implies that 
the boson field is not dynamical at the UV scale Ap^ 
Its dynamics as “measured” by the wavefunction renor¬ 
malization Aip^k is solely generated by quantum fluctua¬ 
tions. Our choice for the initial conditions ensures that 
our ansatz for the effective average action F^, is identical 
to the well-known partially bosonized action at the 
UV scale A, as it should be. 

By applying the Wetterich equation (13) to the 
ansatz (15), the RG flow equations for the various cou¬ 
plings can now be derived. Details as well as our choices 
for the regulator functions can be found in App.f^and[B| 
In the following, we only discuss the general properties 
of these equations. 

The general form of the flow equation for the Yukawa 
coupling reads 


hdkhp — ^^'qA^kkp . (23) 

Note that we have dropped terms on the right-hand side 
which are only non-zero in the regime with broken C/(l) 
symmetry but vanish identically otherwise. Since we are 
primarily aiming at a determination of the phase bound¬ 
aries, but not at a quantitative prediction of a particular 
observable (e.g. the Bertsch parameter) for which these 
terms may become important, we expect this approxi¬ 
mation to be justified. The flow of the Yukawa coupling 
is then purely driven by the anomalous dimension of the 


In our numerical studies, we have chosen a finite value for A = 
lOOOydl 3-ucl therefore a finite value for the initial condi¬ 
tion Ap^k—A = 1 such that it is consistent with Eq. \22\. 


boson field, which be can be split into two distinct con¬ 
tributions: 


lb , ip 

r]A,k = Vam + %.fc • 


(24) 


The contribution rf^ ^ is built up by one-particle irre¬ 
ducible (IPI) diagrams with no internal boson lines but 
only fermion lines. On the other hand, f. receives con¬ 
tributions from IPI diagrams with at least one internal 
boson line and, in our present truncation, it is found to be 
non-zero only in the regime with broken C/(l) symmetry. 

Finally, the flow of the renormalized effective order- 
parameter potential can be conveniently split up into 
three terms: 


kdkUk = VA,kpU'k + [kdkUkf + [kdkUk] 


(25) 


The first term accounts for the renormalization of the 
field (^, the second term is simply a pure fermion loop, 
and the third term is nothing but a pure boson loop. 
The last two are both dressed with suitable regulator 
insertions (see also Appendix 11- 

Using the initial conditions detailed above and inte¬ 
grating the coupled system of flow equations given in 
Eqs. (23)-(25l, we can extract values of physical observ¬ 


ables such as the Bertsch parameter, the fermion gap, 
or the critical temperature. As in the mean-held case, 
a nontrivial global minimum at Uk=oiPo) signals spon¬ 
taneous 17(l)-symmetry breaking associated with a h- 
nite order parameter po = PoPo or, equivalently, a hnite 
fermion gap Ag = hlpo- 

Due to the highly nonlinear, coupled structure of the 
system of how equations, the solutions need to be found 
numerically. As discussed in Sec. |IVB| above, a mean- 
held analysis suggests the existence of hrst-order phase 
transitions. Therefore, we choose to discretize the k- 
dependent ehective potential on a grid in held space (see 
also Ref. [31] for details on our numerical implementa¬ 
tion). 

Before discussing the numerical results from our RG 
how equations, we briehy discuss more precisely in what 
sense Eqs. (23 )-([2^ represent an extension of the conven¬ 
tional mean-held approximation of Sec. EYi Diagrammat- 
ically, Eq. (13) has a simple one-loop structure, where the 


internal lines of Feynman diagrams are given by the so- 
called full propagators ~ (F^^^ -|- Rk)~^- The mean-held 
approximation, on the other hand, takes into account 
Feynman diagrams only with internal fermion lines, but 
with no internal boson lines. In our RG approach, this 
implies that the how equation for the wavefunction renor¬ 
malization and for the effective potential simplify consid¬ 
erably: 


1lA,k — PAk ’ 


(26) 


and 


kdkUk = VA,kpUk + [kdkUk]'^ , 


(27) 
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where, diagrammatically, 


[dkUkf 


s"'. n .■ 




(28) 


n—0 


Moreover, we observe that the flow equation for the ef¬ 
fective potential reduces to 


kdkUki^) = [kdkUkm] 




(29) 


if we rewrite it in terms of the unrenormalized anologue 
of p, namely p = p/A^^k- For the unrenormalized Yukawa 
coupling, we find dkhk = 0. Thus, in this limit the RG 
flow of the effective order-parameter potential Uk and 
wavefunction renormalization are no longer coupled, 
and so the flow of the Yukawa coupling is trivial. This 
implies that the flow equations of Uk and A^^^k can be 
solved independently. Diagrammatically, we find that a 
fermion loop with two external bosonic lines attached to 
it contributes not only to the flow of Uk (see Eq. (28l), 
but also to the flow of A^^k, be. to The con¬ 

sequences of this observation will be discussed in detail 
below. At this point, we note that we do indeed recover 
the well-known mean-field solution for the effective po¬ 


tential at Q = 0 from integrating Eq. (291 with respect 


to the RG scale k, as one should. Here, the restriction 
on the case <5 = 0 comes from the need to project onto 
constant cp in order to obtain explicit expressions for the 
flow equations (24) and (25). Thus, at the level of our 


present approximations, our RG approach is unable to 
resolve inhomogeneous phases explicitly. However, our 
present setup is not “blind” to the fact that inhomoge¬ 
neous condensates may emerge and govern the ground- 
state dynamics, as we discuss next. 


B. Results 


1. Flow of the fermion loop: mean-field and beyond 


As alread y m entioned above, integrating Eq. (28) (or 
rather Eq. (B2)) and minimizing Uk=o{^‘^) yields the 


well-known mean-field phase diagram of mass- and spin- 
imbalanced Fermi gases for the case where the possibility 
of inhomogeneous phases was not taken into account ex¬ 
plicitly (see e.g. Refs. [HlllIlliTll^). However, as men¬ 
tioned above, a fermion loop with two external bosonic 
legs contributes not only to the flow of Uk, but also to 
r]A,k, see Eq. (24). Recall that the flow equations for Uk 
and riA,k in the mean-field approximation are decoupled, 
implying that the fermion gap Aq = h'^po = h^po is in¬ 
dependent of A^p^k■ Although riA,k has therefore no direct 


Note that A^ k directly contributes to the inverse boson propaga¬ 
tor being the second functional derivative of the effective action 
with respect to the boson fields. 



Figure 4. Mean-field phase diagram at T — 0 showing the re¬ 
gion of vanishing (red/dark-shaded area) bounded 

by the green/dotted line. For large m > 0.6, this line associ¬ 
ated with inhomogeneous pairing (see main text) agrees very 
well with the (red) dot-dot-dashed line depicting the lower 
bound for bound-state formation with P > 0 as obtained 
from our two-body analysis, see Sec. m and Fig. The 
orange/dashed line depicts the transition line from the inho¬ 
mogeneous phase to the normal phase as obtained with a FF 
ansatz, see also Figs. and 


impact on the value of Aq, the flow of riA,k still contains 
important information about ground-state properties. 

Indeed, solving the flow equation for alongside 

with the one for Uk, it turns out that the long-range 
limit (fc —)■ 0) cannot always be reached due to a pe¬ 
culiar behavior of the RG flow. This is seen in Fig. 
where we show a red/dark-shaded region bounded by a 
green/dotted line. Within that region, becomes zero 
for a finite value k = fcbreak, (see inset of Fig. for illus¬ 
tration). Strictly speaking, the RG flow breaks down at 
this point. Had we only considered the flow of Uk, we 
would have recovered the well-known mean-field phase 
diagram and easily overlooked this instability, which oc¬ 
curs at the same level of truncation in terms of Feyn¬ 
man diagrams. On the other hand, this instability can 
be expected to be an artifact of our ansatz (15), which 
yields Eq. (20) for the boson propagator. The latter re¬ 
quires Af^^k > 0 in order to be physically meaningful. If 
the boson propagator evaluated at po (or, equivalently, 
Po) turns out to be not positive (semi-)definite, then the 
configuration pq extracted from the computation of Uk 
cannot be the true ground-state of the theory. Of course, 
this instability does not really come as a surprise. In 
Sec. |IV[ we have already shown that, for a simple FF- 
ansatz, a regime governed by an inhomogeneous conden¬ 
sate exists in the phase diagram. Within this regime, 
a homogeneous condensate does not represent the true 
ground state, but rather an “excited” state in the spec¬ 
trum of the theory. 

In our present RG study, we do not employ a spe¬ 
cific ansatz to parametrize the spatial dependence of the 
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Figure 5. Sketches of possible functional forms of the in¬ 
verse bosonic propagator in terms of spatial momentum. The 
green/solid line corresponds to a homogeneous phase, whereas 
the red/dashed and blue/dotted lines sketch two cases with 
a nontrivial minimum at Qmin ~ 0(Q) > 0 favoring the dom¬ 
inance of modes associated with bound states with a finite 
center-of-mass momentum. In case of the red/dashed curve, 
this is linked to a negative 


condensate, such as an FF ansatz. Rather, we study 
the stability of the RG flow under the assumption that 
the condensate does not exhibit spatial dependence. The 
momentum structure in Eq. (20) can be seen as the first 
nontrivial term in a derivative expansion of the exact 
propagator about zero momentum. The fundamental as¬ 
sumption for this ansatz to be valid is that bosonic modes 
with spatial momentum |g| « 0 contribute predominantly 
to the flow. For an inhomogeneous phase this is obvi¬ 
ously not the case, as modes with momenta around \Q\ 
are rather expected to dominate the dynamics. The in¬ 
verse boson propagator close to k — fcbreak is therefore 
expected to be shaped as the (red) dashed line in Fig. 
instead of the green (solid) curve which corresponds to 
Eq. (201 with Ap,^k > 0 (see also Ref. [SS]). It might 
therefore be tempting to associate the existence of a fi¬ 
nite scale fcbreak with the formation of an inhomogeneous 
condensate directly. For the following reasons, however, a 
sign change of should only be considered as a “hint” 
for the emergence of an inhomogeneous condensate: 


• Only for A: —> 0 (long-range limit), strictly speak¬ 
ing, results extracted from the RG flow are phys¬ 
ically observable. Here, we find f 0 at 

fcbreak ^ it may be the case that Ap,^k —f 0 

is negative only for a finite A:-range below fcbreak and 
then becomes positive again in the limit fc —> 0. 
This behavior of Ap,^k may be considered as in¬ 
dicating the existence of an “inhomogeneous pre¬ 
condensation” phenomenon, in analogy to the well- 
known precondensation effect at finite tempera¬ 
ture [Ml and finite spin imbalance |31| . As 
discussed in Ref. [65], it is necessary to employ a 
modified ansatz for the momentum dependence of 
the inverse boson propagator for scales fc < fcbreak 
in order to reach the long-range limit, fc —)■ 0. Oth¬ 


erwise, it is a priori unclear whether we indeed have 
A^^k < 0 for all k < fcbreak, i-e. if the dominance 
of modes favoring an inhomogeneous ground state 
persists in the deep IR. An investigation of this is¬ 
sue is beyond the scope of the present work. 

• The agreement between the boundaries derived 
from our two-body study on the one hand, and the 
criterion associated with the appearance of a zero 
in the flow of on the other, is stunningly good 
for fh > 0.6 (see Fig. |^. It can already be deduced 
from Fig. that inhomogeneous pairing is partic¬ 
ularly preferred for such highly mass-imbalanced 
systems. In fact, our conventional mean-field study 
with a FF ansatz already suggested the existence 
of an inhomogeneous phase in this regime (see the 
domain enclosed by the orange/dashed line and the 
black line in Fig. |^. However, we also observe 
that bound states with finite center-of-mass mo¬ 
menta as well as the regime with fcbreak > 0 are 
found to exist well beyond this FF-type phase. This 
discrepancy may be traced back to the fact that 
a simple FF ansatz may be insufficient in some 
parts of the phase diagram. However, this can also 
be viewed as a manifestation of the existence of 
bosonic bound states (or the dominance of finite- 
momentum bosonic modes), which does not neces¬ 
sarily imply condensate formation, even at T = 0. 

• For h < 0, the domain bounded by the crite¬ 
rion of A^^k 0 extends down to m « 0.22 for 
h = — 1 , whereas no bound states with finite center- 
of-mass momentum are found below rh = 0.46 and 
h = —0.86, see Figs. and The discrepancy be¬ 
tween these two lines should not be too surprising 
for small h. In fact, the behavior of Ap,^k is only 
indirectly linked to the appearance of bound states. 
It is reasonable to expect that the momentum de¬ 
pendence, i.e. the functional form, of the propaga¬ 
tor is dominated by the existence of deeply bound 
two-body states as they are present for fh > 0.6, 
where the two lines are in very good agreement. In 
this regime, condensation of these states may in¬ 
deed be energetically favored and the emergence of 
an instability in the RG flow associated with the 
observation Ap,^k 0 may then be viewed as a 
strong indication that the formation of an inho¬ 
mogeneous condensate is most favorable. However, 
once the binding energy becomes smaller, the func¬ 
tional form of the propagator becomes progressively 
more dominated by modes with momenta signif¬ 
icantly different from the center-of-mass momen¬ 
tum of the lowest lying two-body bound state for 
given values of h and fh. Even if there is no bound 
state at all in our analysis of the two-body problem, 
there is no reason why bosonic fluctuations in gen¬ 
eral, and those with finite momentum in particu¬ 
lar, should become irrelevant. Consequently, in do¬ 
mains with only shallowly bound two-body states. 
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as it is the case for decreasing h, we may expect at 
best qualitative agreement between our predictions 
from the two-body problem and the study of the 
RG flow of 

In summary, —>■ 0 signals the dominance of finite- 

momentum bosonic fluctuations at least for a certain 
range of values of the RG scale k < fcbreak- Strictly 
speaking, however, even positive values of A^p^k^o do not 
imply that a low-order derivative expansion of the effec¬ 
tive action captures the correct dynamics with respect to 
inhomogeneous phases. In fact, the true inverse propaga¬ 
tor may, e.g., be shaped as depicted by the blue/dotted 
line in Fig. In any case, it is also clear that it is not 
possible to discriminate irrevocably between normal, ho¬ 
mogeneous and inhomogeneous phases only with the aid 
of the criterion Ap,^k —t 0. If, on the other hand, an 
existing homogeneous condensate (as found in a study 
of the RG flow of the potential Uk) overlaps with a do¬ 
main where A^p^k —>■ 0 is also observed in the flow, this 
is a strong indication of an inhomogeneous phase in that 
same domain. This point will be detailed further in the 
next section, where we present results from the RG flow 
including bosonic fluctuations. 


2. Full Flow Equations: Impact of Bosonic Fluctuations 


Let us now study the effect of the bosonic contributions 
r]'^ j, and [dkUk]‘^ in our set of flow equations. Due to the 
derivatives oi Uk with respect to p in the flow equations 
for Ap^k and Uk, and the nontrivial dependence of these 
equations on hp (see Eqs. (B3) and (Bll)), the flow equa¬ 


tions (23)-(251 are now coupled in a highly nontrivial way. 
As discussed in Sec. |V A[ we solve this set of flow equa¬ 
tions numerically. However, it does not appear feasible to 
either start the integration at fc = A —> oo, nor to reach 
the IR limit /c —>■ 0 exactly. In the IR regime, we have 
stopped the RG flow at fciR « lQ~^y/jL in those domains 
of parameter space that are close to a second-order phase 
transition. The remaining uncertainty in the critical tem¬ 
perature, as measured by the fermion gap A^, turned out 
to be well below the percent level. In domains close to 
a first-order transition, the introduction of a finite IR 
cutoff fciR is even less severe as the order parameter 
is found to jump typically on scales k ^ O{,/Jl/10) (see 
also Ref. [3T] for a detailed discussion of the numerical 
implementation in the mass-balanced limit). 

The phase diagram from our RG study including order- 
parameter fluctuations is shown in Fig. We first note 
that the qualitative features of the phase diagram are 
similar to the mean-field phase diagram discussed above. 
However, quantitative corrections are found to be sizable. 
For rh = h = 0, for example, we find that the critical tem¬ 
perature T is significantly changed from « 0.66 in the 
mean-field approximation to « 0.40, in good agree¬ 
ment with recent experiments [68) . Monte Carlo stud¬ 
ies as well as with recent RG studies [SB Eg. 



Figure 6. Phase diagram obtained from the full set of 
flow equations (|23|l-( 25 I including bosonic fluctuations. The 


red/dark-shaded area marks the domain where an inhomo¬ 
geneous phase is likely to exist, see main text for a detailed 
discussion. 


Lowering h, starting from h = ffi for fixed m < 0.53, 
the order of the phase transition changes from second 
to first along a critical line Tc{h,m). In this ffi-regime, 
we find Tcp(h,rh) « 0.19...0.20 for the critical point 
where the nature of the transition changes from second to 
first order. For temperatures below T^p{h,m) and given 


m, the superfluid phase is enlarged in the /i-direction 
compared to the mean-field result, as also found in the 
mass-balanced case [3T]. For T <C Tc the integration of 
Eqs. (23l-(|2g) becomes numerically more and more in¬ 


volved, as the right-ha nd s ides b ecom e discontinuous at 
T = 0 (see also Eqs. ( |B2[ ) and (B6)). Contrary to the 
mean-field case, we therefore do not present results for 
the strict zero-temperature limit here but restrict our¬ 
selves to temperatures T > 0.1. However, we expect 
the position of the phase transition lines in the zero- 
temperature limit to still be in reasonable agreement with 
the ones obtained for T = 0.1. 


Lowering h, starting from h = m for fixed m > 0.53 
and T < Tjfp(m) with Tjp(m) « 0.19...0.21, we find 
that a critical value exists at which A^^k 

tends to zero, potentially indicating a transition from 
a homogeneous superfluid phase to an inhomogeneous 
phase. In fact, we find that Ap,^k tends to zero at a finite 
RG scale fcbreak for all h G [—1,/ibreak(di, T)]. As dis¬ 
cussed in Sec. |VB 1[ this does not necessarily imply that 
there exists an inhomogeneous phase for the whole do¬ 
main h G [—1, /ibreak(fh, T)] for a given to and T. There¬ 
fore the red/dark-shaded area in Eig. [^only represents a 
domain where an inhomogeneous phase is most likely to 
be found. Our mean-field study supports this interpreta¬ 
tion: In Fig. we find reasonable agreement for large to 
between the onset of inhomogeneous condensation and 
the line defined by the largest value of h for which A,p^k 
still tends to zero at a finite scale /cbreak, for a given value 
of TO. For this reason, we consider our A,^_fc-criterion to 
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provide a reasonable estimate for the boundary separat¬ 
ing the phase with a homogeneous condensate from the 
inhomogeneous phase, even beyond the mean-field limit. 
However, within our present setting, it is not possible to 
detect the transition from the inhomogeneous phase to 
the normal phase. The associated phase-transition line 
is therefore not shown in Fig. Instead, we indicate 
the presumed existence of this phase-transition line by a 
fading out of the red/dark-shaded domain p] 

In our mean-field studies, conventional and RG, we 
have considered the zero-temperature limit explicitly. 
This allowed us to compare directly the lower bound 
for pairing with finite center-of-mass momentum and the 
bound from our H,^_fc-criterion. Since the strict zero- 
temperature limit is difficult to reach from a numeri¬ 
cal point of view when bosonic fluctuations are taken 
into account, we have restricted ourselves to tempera¬ 
tures T > 0.1. However, this renders a direct comparison 
with the results from our two-body analysis impossible. 
Moreover, taking into account bosonic fluctuations, our 
^(, 5 ,fe-criterion used to detect the onset of inhomogeneous 
phases has to be taken with some care: In Fig. the 
fermion gap Aq is depicted as a function of h for two 
different temperatures and fixed rfi = 0.74 correspond¬ 
ing to a mixture of ®Li and For f = 0.225, the 
green/dot-dashed line exhibits the typical behavior ex¬ 
pected for a second order transition to the normal phase. 
Only slightly below this temperature, at T = 0.20, the 
behavior of the condensate is now found to be quite dif¬ 
ferent. Already before reaching the regime defined by 
h < /ibreak, where tends to zero in the RG flow 
(indicated by the red/dark-shaded area), the gap Aq is 
found to increase rather than decrease. Such a behavior 
is not seen in the mean-field approximation. In fact, for 
any finite T and fixed rh, we find that Aq is a strictly 
concave function. One may speculate that the observed 
increase in the gap is due, e.g., to numerical artifacts. 
Indeed, the value of Ag is more sensitive to numerical in¬ 
accuracies in this region than anywhere else in the phase 
diagram and should therefore be taken with some care. 
However, we have checked that the general observation 
of an increasing gap cannot be traced back to numerical 
instabilities. 

The inset of Fig. hints at the true reason for the 
unexpected behavior of the gap as a function of h for 
given m at low temperatures: It shows the fc-evolution 
of (red/solid line) and A^ (blue/dotted line) for a 
point in a region of the phase diagram where the gap 
increases again but A^^^ is still positive on all scales k. 
In this region (as exemplified by the inset), we observe 
that Aip^k undergoes a sharp decrease. From the depen¬ 
dence of the RG flow of the Yukawa coupling on A^^k 


Of course, a priori, it is unclear whether a transition from an 
inhomogeneous to a normal phase exists after all. However, our 
mean-field study in Sec. 113 our analysis of the two-body problem 
in Sec. |III[ and studies of the Fermi polaron cziiaffl suggest 
that such a transition indeed exists. 



Figure 7. Fermion gap for two different values of T and fixed 
m = 0.74 (corresponding to a Li-K mixture) as a function 
of h. The red/dark-shaded area indicates the regime where 
A,p,k —t 0 for T = 0.20. The inset shows the fc-evolution of 
and for {T,h) = (0.20,0.51). 


via the anomalous dimension 77 ^ ^ (see Eq. (23)), it is 
clear that h^p experiences a sharp increase when A^^k 
decreases. The sharp increase of the Yukawa coupling 
also increases the gap up to values larger than those in 
which no (strong) decrease of A^^k is present at all (e.g. 
deep in the homogeneous superfluid phase). In the inset 
of Fig. we also observe that the increase of the gap 
induced by the decrease of A^^k can potentially counter¬ 
balance the decrease of A,p^k- This effect is intimately 
related to the presence of bosonic fluctuations. In fact, 
the strong increase of the gap leads to a suppression of 
Feynman diagrams with internal fermion lines. In par¬ 
ticular, purely fermionic diagrams present in mean-field 
studies are parametrically suppressed by a (strong) in¬ 
crease of the gap. Thus, once a gap is generated, the 
RG flows of the potential, the Yukawa coupling, and the 
wavefunction renormalization A^^k are mainly driven by 
purely bosonic diagrams (without any internal fermion 
lines). As (purely) bosonic and fermionic contributions 
come in general with opposite signs, bosonic fluctuations 
tend to counterbalance the decrease of A,p^k as well as the 
increase of the gap, which is induced by purely fermionic 
diagrams (see also Appendix for the RG flow equa¬ 
tions of the various quantities). Since the class of purely 
fermionic diagrams is the only one present in the mean- 
field approximation, this explains why such a counterbal¬ 
ancing effect is not observed in our mean-field study. 

The observed decrease in A,p^k can at least partially 
be traced back to the fact that a change of the spin im¬ 
balance parameter h induces a mismatch in the Fermi 
momenta and Fermi energies associated with the spin- 
up and spin-down fermions. Within our RG framework 
this implies (loosely speaking) that the two spin com¬ 
ponents contribute to the RG flow over different ranges 
of scales k. In any case, the decrease in actu¬ 

ally becomes stronger when we decrease h for a given fh 
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starting from the line of equal Fermi momenta h = m. 
Eventually, we observe that the mismatch between the 
Fermi momenta becomes so large that the associated 
decrease of can no longer be compensated by the 
presence of bosonic fluctuations, resulting in t 0 

for h < hbreak(fh, T) at a given finite RG scale fcbreak, at 
least within our present truncation. 

We close this section with a word of caution regard¬ 
ing the increase of the gap towards smaller values of h 
in some domains of the phase diagram (see Fig. [^. At 
present, it is unclear whether this behavior is a phys¬ 
ical effect and traces of it should be expected in fu¬ 
ture experimental studies, or whether it is “cured” when 
considering a suitable extension of our present trunca¬ 
tion. For example, a resolution of the full momentum 
dependence of the inverse boson propagator might 
be required to resolve this issue since, in the relevant 
domain, may assume a form as exemplified by the 
blue/dotted line in Fig. However, a detailed analysis 
of this issue is left to future work. Here, we restrict our¬ 
selves to conclude that our RG study already suggests 
that bosonic fluctuation effects are of great importance 
in studies of mass- and spin-imbalanced unitary Fermi 
gases. We have found that the extent of the various 
phases in parameter space (normal, homogeneous super¬ 
fluid, inhomogeneous superfluid) can change significantly 
relative to mean-field studies. In particular, our RG anal¬ 
ysis indicates that the size of the inhomogeneous phase 
is extended to smaller values of m. More specifically, our 
mean-held studies suggest that m > 0.65 (m^/m^ > 4.7) 
is required to form an inhomogeneous condensate. Re¬ 
call that fh ~ 0.74 (m^/m^ > 6.7) for a Li-K mixture. 
Taking bosonic fluctuations into account, we hnd that in¬ 
homogeneous condensates may already appear for mass 
imbalances fh > 0.53 > 3.2). 


VI. SUMMARY 

We have studied the phase diagram of mass- and spin- 
imbalanced unitary Fermi gases in three dimensions, with 
an emphasis on the detection of inhomogeneous conden¬ 
sates. To this end, we have employed various approaches. 
Since the formation of two-body bound states can be 
viewed as a necessary condition for condensation, we have 
analyzed the two-body problem in the presence of (inert) 
Fermi spheres associated with the two spin components, 
and computed the binding energy as well as the center-of- 
mass momentum of the lowest lying bound state. This 
helped us guide and analyze our studies of the many- 
body phase diagram. Indeed, we have found that our 
two-body bound-state problem allows us to understand 
many features of the many-body phase diagram, such as 
the emergence of a domain in parameter space governed 
by a homogeneous condensate. Moreover, our study of 
two-body bound states suggests the existence of a region 
in the many-body phase diagram where an inhomoge¬ 
neous condensate is formed, see Fig. 


Our explicit studies of the many-body problem indeed 
confirm the qualitative picture of the phase diagram sug¬ 
gested by the two-body problem. In particular, we have 
only found indications of inhomogeneous phases in those 
regimes of the many-body phase diagram in which also 
the center-of-mass momentum of the lowest-lying two- 
body bound state is Hnite. The agreement of the phase 
boundary of the homogeneous superfluid and inhomoge¬ 
neous phases with the lower bound for pairing with finite 
center-of-mass momentum is only qualitative when we 
consider a Fulde-Ferrell ansatz for the inhomogeneity in 
our mean-field study. On the other hand, the agreement 
of this lower bound with the instability linj^ associated 
with the homogeneous-inhomogeneous transition is stun¬ 
ning for large fh > 0.5 in our RG mean-field study, which 
does not rely on an explicit ansatz for the inhomogeneity. 
Taking bosonic fluctuations into account in our RG study, 
the latter two lines are still in accordanceE3 For mass im¬ 
balances fh < 0.5, the analysis of the two-body problem 
still suggests the existence of a regime with bound states 
with finite center-of-mass momentum. Our many-body 
studies, however, suggest that inhomogeneous conden¬ 
sates are unlikely to be found in this regime. This can be 
traced back to the fact that only shallowly bound two- 
body states exist there. Moreover, it also indicates that 
macroscopic condensation of such bound states is not en¬ 
ergetically favored. 

For large fh > 0.5, on the other hand, a regime with 
deeply bound two-body states with finite center-of-mass 
momentum can be identihed, suggesting that the forma¬ 
tion of an inhomogeneous condensate is energetically fa¬ 
vored. Indeed, we find that this regime overlaps with 
the regime where homogeneous condensation is found in 
our many-body analysis whenever inhomogeneous pair¬ 
ing is not taken into account. Therefore, the instability 
line found in our RG study including bosonic fluctuations 
may indeed be identihed with the transition line from a 
superhuid phase with a homogeneous condensate to one 
with an inhomogeneous condensate, at least for fh > 0.5. 
Finally, for small mass imbalances fh ^ 0.25 we do not 
hnd any indication of inhomogeneous phases, neither in 
our analysis of the two-body problem nor in our many- 
body studies. 

In summary, our RG study beyond the mean-held limit 
suggests signihcant modihcations to the mean-held ex¬ 
tent of the various phases in parameter space. For ex¬ 
ample, for the mass- and spin-balanced Fermi gas, we 
hnd that the critical temperature Tc agrees well with 
state-of-the-art Monte Garlo studies, as already found 


Recall that, for a fixed temperature, this line is defined by the 
largest value of h for which still tends to zero at a finite 

scale fcbreak for S' given value of fh. 

Note that a direct comparison of both lines is not possible at this 
stage since we restricted ourselves to temperature T > 0.1 in our 
RG study including bosonic fluctuations. As discussed above, a 
study of even lower temperatures is numerically challenging and 
costly and is therefore left to future work. 
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in recent RG studies [El m\- For small mass imbal¬ 
ances, moreover, we find that the phase with a homo¬ 
geneous condensate is extended to larger spin imbal¬ 
ances \h\ in agreement with a previous RG study of 
the mass-balanced case m- For large mass imbalances, 
however, we observe that the extent of the phase with 
a homogeneous condensate shrinks in the h direction. 
Moreover, we find strong indications for the emergence 
of an inhomogeneous phase in this regime. The ex¬ 
tent of the latter phase in the /i-direction actually in¬ 
creases for increasing mass imbalance. Our estimate for 
a lower bound on m for the occurrence of an inhomo¬ 
geneous condensate is m « 0.53. This value of m cor¬ 
responds to a mass ratio m^/m^ « 3.2 which is con¬ 
siderably lower than the mass ratio associated with a 
Li-K mixture « 6.7). Since three-body effects 

are expected to be significant for mass ratios associated 
with the Li-K mixture and beyond (see, e.g.. Refs. [69l- 
EI]), future experimental searches for, e.g., signals of in¬ 
homogeneous condensation in the regime m^/m^ > 6.7 
are highly challenging. In this respect, our estimate for 
the lower m-bound for the occurrence of inhomogeneous 
phases may help to guide future experiments. Our results 
suggest the existence of a much larger regime in parame¬ 
ter space where inhomogeneous condensation may be de¬ 
tected with a new variety of mixtures of fermion species 
other than Li-K, but without suffering from three-body 
effects. 


Note that in both cases only spatial momenta are regu¬ 
larized which allows us to perform the Matusbara sums 
in the loop integrals analytically. 

The bosonic regulator function is the same as in 
previous studies [iS]. Only the argument has been mod¬ 
ified by a factor (1 — ?7i^) to accommodate for the gener¬ 
alized dispersion relation in the mass-imbalanced case. 

Intuitively, the two fermion species should be regulated 
separately, since their dispersion relations differ by con¬ 
struction in the imbalanced case. This appears to be 
particularly important since the presence of Fermi seas 
in general requires fluctuations below and above the re¬ 
spective Fermi levels to be treated differently. However, 
the analytic computation of the frequency sum reveals 
that in the critical cases, only an “averaged” Fermion 
dispersion appears in the RG flow equations. Thus, the 
“averaged” regulator (AI) is sufficient to render the in¬ 
tegration over spatial momenta finite. The structure of 
the flow equations is thereby greatly simplified. In con¬ 
trast to the mass-balanced but spin-imbalanced case m, 
to our knowledge, an analytical computation of the mo¬ 
mentum integrals is not possible. We add that IR observ¬ 
ables from a RG flow controlled by the regulator func¬ 
tion (Al) on the one hand and by separately regulated 
fermion propagators on the other hand has found to agree 
extremely well in the mass-balanced case, see Ref. m- 
We therefore expect that our choice (Al) can safely be 
employed in the mass-imbalanced case as well. 
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In this appendix we provide explicit expressions for the 
flow equations of the effective potential and the boson 
anomalous dimension rjA,k- Since they can be derived 
largely along the lines of the balanced or spin-imbalanced 
cases, we focus mainly on the peculiarities of the mass- 
imbalanced case. For further details, we refer the reader 
to earlier work, see, e.g., Refs. [Ell45ll48] . 


Appendix A: Regulator Functions 


Apart from the very defining properties of a regulator 
function Rk in the fRG formalism (see Ref. |6^), there 
is considerable freedom in the choice concerning the pre¬ 
cise shape of Rk which even allows to optimize the RG 
flows, see, e.g., Refs. EaES]. Guided by these powerful 
principles of optimization, we employ the following set of 
functions for the fermion and boson fields, respectively: 

Rf (z) = (sign(z) -z)0{l-\z\), z= ^ ’ (^1) 


and 


Kiv) = - y), y 



(A2) 


1. Flow of the effective potential 

The flow equation of the order-parameter potential 
Uk{p) is obtained from 

Note that, for convenience, we have written the com¬ 
plex boson field as a sum of its real and imaginary 
parts, if{T,x) = Moreover, we 

add that, due to the the fact that Uk only depends 
on p = ip*ip, the projection on ipi = 2y/o and p 2 = 0 
can be applied without loss of generality. Eventually, 
this yields 

dkUk = \-nA,kpU'k + [dkUkf + [dkUkY (Bl) 

with 
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[dkUk]'^ = ^ [ dz ^ V aNp (m{z + fi)-h + a^/l + ws) , 

27’" Jmax[-/i-l] V / 


(B2) 


[afeC/fe]‘" = 


V2k^ 


Stt^ (1 — m^) ^ 


i'-f) 


1 + wi 11 

1 + W2 


■ W2 


1 + Wl 


- + A^g (yi + wi\/l + ^ 2 ) 


(B3) 


Here, A^p and A^b correspond to the Fermi and Bose- 
Einstein distribution functions, respectively: 


r 


the “steepness” of the Fermi function makes the numer¬ 
ical treatment increasingly challenging as T is lowered. 


Np{x) = 


1 

^/T + 1 ' 


Np, = 


1 

- 1 


(B4) 


Furthermore, the following quantities have been intro¬ 
duced: 


Uk 


U’k + '^pUk 


^3 = ^ ■ (B5) 


Quantities Q divided by a factor are written as Q, e.g. 
jl = . 


Since Aj, = Ktp = k'^w^, Eq. (B2) does not depend 


solely on p. This allows us to identify the solution from 
standard mean field theory with the solution of the purely 
fermionic flow of in the limit k —>■ 0, as discussed at 
the end of sect. |V A[ Furthermore, we note that 


lim Np(x) = 0(—x). 

T-J-O 


(B6) 


Thus, the right-hand side of Eq. (B2) becomes non- 


analytic in the zero-temperature limit. While the flow 
of C/fe in principle exists for arbitrarily small but finite T, 


2. Boson anomalous dimension 

In order to extract the running of the wavefunction 
renormalization parameter A^p^k via the boson anomalous 
dimension r]A,k = —kdklaAp,^k, the following projection 
rule is employed: 


77A,fc = 


d 


Ap^k l-m? dq"^ 


9k [(-fV)i2(0’^] 


J9=0 


where 


_ •0,1 I 0,2 I ip 

= ^A,k + VA,k + klA,k > 


{Pv)i2 


(B7) 


TkT^ 


Sipi{-p) Sip2iq) 


(B8) 


p*=-p=^2=0,ipi = y/2poA 


Note that we do not consider a p-dependent pA,k- In¬ 
stead, riA,k is projected onto the /c-dependent minimum 
po,fc of the effective potential. The explicit expressions 
for the different contributions to r]A,k are given by: 


0,1 


— 2 ^ _ 




{fl +k)^ 9{fl +k) 


m 67r2/c(l-pu;3)2 

aNp {m{p + k) — h — a^l -|- -|- + w^Np {m{p + n) — h — a^/l + 


(B9) 


j/j.2 

PA,k = 


dz 


{z + p)i 


1^3(7Ap (^m{z + p) — h + tr-v/T+dca^ 


67r2A:(l - m2) (l-Pwaji Y±\ 

-3-\/l -I- W'iNp (jh{z + p) — h — a^AA^A^I}^ — a{l + w^)Np {m{z + p) — h — a^Jl + 


(BIO) 


^ PoC/f V2 
Vlk = 


(1 — m2)2 37r2fc -p ra^))! -p ^ 2 )] 2 


2 [1 -k 2A7 b {VI + WWI + W 2 ) - 2A(5 (v'l + wiVl + W2)] . (Bll) 


Here, the primes denote derivatives of the thermal dis- Np^^{x) = dNp/^{x)/dx, and the Wi's are understood 
tribution functions with respect to their arguments. 
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to be evaluated at p = po,fc- 

The contributions in Eqs. (B9) and (Bll) represent 
straightforward extensions of their mass-balanced ana¬ 
logues. The situation is different for 77 ^’^ since it is not 
even present in the mass-balanced case. It becomes im¬ 
portant only at relatively large mass imbalances since it 
is proportional to fh? / (1 — fh^). 

Taking a closer look at the bosonic contribution 77 ^ 
it becomes clear why a finite po,fc may have such an enor¬ 


mous influence on the RG flow of the theory (see also 
Fig. 0 : k vanishes identically in the 17(l)-symmetric 

regime, i.e. for po,fe = 0 - 

It is not too surprising that bosonic fluctuations are 
generally amplified by mass imbalance. This becomes 
apparent by the fact that both 77 ^ and [dkUkY are pro¬ 
portional to (1 — fh?)~i. A detailed analysis how this 
affects the properties and the extent of inhomogeneous 
phases for large fh is left to future work. 
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